function [output, setup] = GPOPSRobustness(varargin)

    addpath('../Densities')
    parser = inputParser;
    parser.addParameter('version',[]) % optimum,
                                      % default_capital_tax, equipment_equal_robot_tax
    parser.addParameter('calibration_version',[])
    parser.addParameter('kappa_R_increase',true)
    parser.addParameter('num_points',50)
    parser.parse(varargin{:})
    auxdata.settings.version = parser.Results.version
    calibration_version = parser.Results.calibration_version

    % whether parameter kappa_R should be increased or decreased, starting
    % from the calibration
    kappa_R_increase = parser.Results.kappa_R_increase;
    
    % number of points for robustness computation
    num_points = parser.Results.num_points;
    
    if kappa_R_increase
        suffix = 'kappa_R_incr';
    else
        suffix = 'kappa_R_decr';
    end
    
    dir_out = '../../../data/generated/matlab/OptimTax/gpops/';
    file_vars = [dir_out,auxdata.settings.version,'_vars_',suffix,'.csv'];
    file_vectors = [dir_out,auxdata.settings.version,'_vectors_',suffix,'.csv'];

    %% Options for debugging
    auxdata.settings.counter = 0; %1, count iterations inside continous
    if auxdata.settings.counter
        counter = 0;
        save('counter.mat','counter');
    end

    %% settings used by GPOPSFunctions and GPOPS preparation
    auxdata.settings.log_w = 1; %1 if w should be transformed to log
    auxdata.settings.standardize = 1; %1 w should be standardized based

    %% initialize gpops_fun object

    par_prod_fun = tools.csv2struct(['../../../data/generated/matlab/Calibration/',calibration_version,'_KrusellRobotsEach_sol.csv']);
    mom_prod_fun = tools.csv2struct(['../../../data/generated/matlab/Calibration/',calibration_version,'_KrusellRobotsEach_mom.csv']);

    par_dist = tools.csv2struct(['../../../data/generated/matlab/Calibration/',calibration_version,'_abilityDistThreeFixC_sol.csv']);
    mom_dist = tools.csv2struct(['../../../data/generated/matlab/Calibration/',calibration_version,'_abilityDistThreeFixC_mom.csv']);

    inequ_par = tools.csv2struct('../../../data/generated/matlab/Calibration/inequ_par.csv');
    constants = tools.csv2struct('../../../data/raw/constants.csv');

    prod_fun = ProductionFunction(par_prod_fun);
    skill_dist = SkillDistribution(par_dist);
    part_dist = ParticipationCostDistribution(par_dist);

    epsilon = constants.epsilon;
    xi = par_dist.xi;
    util_fun = UtilityFunction(epsilon, xi);

    par.r = inequ_par.inequ_par;
    par.revenue_requirement = mom_prod_fun.revenue_requirement;

    par.q_B = mom_prod_fun.q_B;
    par.q_E = mom_prod_fun.q_E;
    par.q_S = mom_prod_fun.q_S;
    par.num_nodes = 100; % can be low if fixing participation
    par.phi_lb = constants.phi_lb;

    auxdata.gpops_fun = GPOPSFunctions(prod_fun,...
                                       util_fun,...
                                       skill_dist,...
                                       part_dist,...
                                       par);

    auxdata.settings.tau_B = par_prod_fun.tau_B;
    auxdata.settings.tau_E = par_prod_fun.tau_E;
    auxdata.settings.tau_S = par_prod_fun.tau_S;

    %% GPOPS preparation
    % bounds for wages (in levels)

    diary
    diary(['gpops_output_', suffix, '.txt'])
    diary on
    wf_levels = constants.wf_gpops;
    for i=constants.w0_gpops
        
        w0_levels = i;
        fprintf('\n=============================\n')
        fprintf('w0_levels = %f, wf_levels = %f\n', w0_levels, wf_levels)
        fprintf('=============================\n')
        
        
        %% approximate density, used for standardization
        if auxdata.settings.standardize
            w_points = linspace(1e-6,500,500);
            Y_M_0 = par_dist.Y_M_0;
            Y_R_0 = par_dist.Y_R_0;
            Y_C_0 = par_dist.Y_C_0;
            
            f_w = auxdata.gpops_fun.skill_dist.f(w_points,Y_M_0,Y_R_0,Y_C_0);
            fit = fit_plognpdf(w_points,f_w);
            auxdata.alpha_fit = fit.alpha;
            auxdata.mu_fit = fit.mu;
            auxdata.sigma_fit = fit.sigma;
        end
        
        if auxdata.settings.log_w==1
            w0 = log(w0_levels);
            wf = log(wf_levels);
            if auxdata.settings.standardize==1
                w0 = (w0 - fit.mu)/fit.sigma;
                wf = (wf - fit.mu)/fit.sigma;
            end
        else
            w0 = w0_levels;
            wf = wf_levels;
        end
        
        theta_M_lb = w0_levels/par_dist.Y_M_0;
        theta_C_ub = wf_levels/par_dist.Y_C_0;
        
        auxdata.w0 = w0_levels;
        auxdata.wf = wf_levels;
        
        % initial values for parameters
        par0.L_M = par_prod_fun.L_M_0;
        par0.L_R = par_prod_fun.L_R_0;
        par0.L_C = par_prod_fun.L_C_0;
        par0.K_B_M = par_prod_fun.K_B_M_0;
        par0.K_B_R = par_prod_fun.K_B_R_0;
        par0.K_B_C = par_prod_fun.K_B_C_0;
        par0.K_E = par_prod_fun.K_E_0;
        par0.K_S = par_prod_fun.K_S_0;
        par0.transfer = par_dist.transfer;
        
        % initial values for control and state
        
        l0 = (par_dist.xi*w0_levels)^epsilon;
        lf = (par_dist.xi*wf_levels)^epsilon;
        c0 = l0*w0_levels; % only used to compute u guesses
        cf = lf*wf_levels; % only used to compute u guesses
        u0 = c0-1/par_dist.xi.*l0^(1+1/epsilon)/(1+1/epsilon);
        uf = cf-1/par_dist.xi.*lf^(1+1/epsilon)/(1+1/epsilon);
        
        % by which factor to scale down
        auxdata.settings.inpscale.l = 1;
        auxdata.settings.inpscale.u = 1;
        auxdata.settings.inpscale.welfare = 1;
        
        l0 = l0/auxdata.settings.inpscale.l;
        lf = lf/auxdata.settings.inpscale.l;
        u0 = u0/auxdata.settings.inpscale.u;
        uf = uf/auxdata.settings.inpscale.u;
        
        %% Initial values, final values and bounds
        smallNum = 1e-3;
        bigNum = 1e3;
        
        % Controls
        lMin                     =      1e-3;
        lMax = (par_dist.xi*wf_levels*2)^epsilon; %allowing for 100% negative
                                                  % top
                                                  % rate
                                                  % States
        uMin                     =      1e-3;
        uMax = lMax*wf_levels - 1/par_dist.xi.*lMax^(1+1/epsilon)/(1+1/epsilon);
        
        
        %-------------------------------------------------------------------------%
        %----------------------- Setup for Problem Bounds ------------------------%
        %-------------------------------------------------------------------------%
        bounds.phase.initialtime.lower      =      w0;
        bounds.phase.initialtime.upper      =      w0;
        
        bounds.phase.finaltime.lower        =      wf;
        bounds.phase.finaltime.upper        =      wf;
        
        bounds.phase.initialstate.lower     =      uMin;
        bounds.phase.initialstate.upper     =      uMax;
        
        bounds.phase.state.lower            =      uMin;
        bounds.phase.state.upper            =      uMax;
        
        bounds.phase.finalstate.lower       =      uMin;
        bounds.phase.finalstate.upper       =      uMax;
        
        bounds.phase.control.lower          =      lMin;
        bounds.phase.control.upper          =      lMax;
        
        results_vars = table();
        results_vectors = table();
        
        par0.L_M = par_prod_fun.L_M_0;
        par0.L_R = par_prod_fun.L_R_0;
        par0.L_C = par_prod_fun.L_C_0;
        par0.K_B_M = par_prod_fun.K_B_M_0;
        par0.K_B_R = par_prod_fun.K_B_R_0;
        par0.K_B_C = par_prod_fun.K_B_C_0;
        par0.K_E = par_prod_fun.K_E_0;
        par0.K_S = par_prod_fun.K_S_0;
        par0.transfer = par_dist.transfer;
        
        
        % compute mass for density normalization
        f = @(w) auxdata.gpops_fun.skill_dist.f(w, par_dist.Y_M_0, par_dist.Y_R_0, par_dist.Y_C_0);
        mass = integral(f,w0_levels,wf_levels);
        auxdata.gpops_fun.skill_dist.mass = mass;
        
        f = @(w) auxdata.gpops_fun.skill_dist.f(w, par_dist.Y_M_0, par_dist.Y_R_0, par_dist.Y_C_0);
        fprintf('\nf(w0) = %f, f(wf) = %f\n',f(w0_levels), f(wf_levels))
        
        % if bounds for parameters are binding in solution, extend
        % bounds and run again
        
        auxdata.settings.consistency_conditions = true;
        auxdata.settings.no_arbitrage_robots = true;
        
        switch auxdata.settings.version
          case 'optimum'
            auxdata.settings.default_capital_tax = false;
            auxdata.settings.default_robot_tax = false;
            auxdata.settings.equipment_equal_robot_tax = false;
          case 'default_capital_tax'
            auxdata.settings.default_capital_tax = true;
            auxdata.settings.default_robot_tax = true;
            auxdata.settings.equipment_equal_robot_tax = false;
          case 'equipment_equal_robot_tax'
            auxdata.settings.default_capital_tax = false;
            auxdata.settings.default_robot_tax = false;
            auxdata.settings.equipment_equal_robot_tax = true;
          case 'sensitivity'
            auxdata.settings.default_capital_tax = false;
            auxdata.settings.default_robot_tax = false;
            auxdata.settings.equipment_equal_robot_tax = false;
        end
        
        auxdata.settings.inpscale.L_M = par0.L_M;
        auxdata.settings.inpscale.L_R = par0.L_R;
        auxdata.settings.inpscale.L_C = par0.L_C;
        auxdata.settings.inpscale.K_B_M = par0.K_B_M;
        auxdata.settings.inpscale.K_B_R = par0.K_B_R;
        auxdata.settings.inpscale.K_B_C = par0.K_B_C;
        auxdata.settings.inpscale.K_E = par0.K_E;
        auxdata.settings.inpscale.K_S = par0.K_S;
        auxdata.settings.inpscale.transfer = par0.transfer;
        
        
        auxdata.settings.fix_participation = false;
        % auxdata.settings.pi_M = 0.722294549151161;
        % auxdata.settings.pi_R = 0.82367788344373;
        % auxdata.settings.pi_C = 0.851776748316647;
        
        auxdata.settings.pi_M = 0.6907;
        auxdata.settings.pi_R = 0.7564;
        auxdata.settings.pi_C = 0.8497;
        
        
        bounds.phase.integral.lower = [...
            -bigNum,...
            -bigNum, ...
            -bigNum,...
            -bigNum,...
            0,...
            1e-6,...
            1e-6,...
            1e-6,...
            smallNum, ...
            smallNum,...
            smallNum];
        
        bounds.phase.integral.upper = [...
            bigNum,...
            bigNum,...
            bigNum, ...
            bigNum,...
            0,...
            1e3,...
            1e3,...
            1e3,...
            bigNum,...
            bigNum,...
            bigNum];
        
        
        bounds.eventgroup(1).lower = [0, 0];
        bounds.eventgroup(1).upper = [0, 0];
        
        bounds.eventgroup(2).lower = [0, 0, 0];
        bounds.eventgroup(2).upper = [0, 0, 0];
        
        if auxdata.settings.default_robot_tax
            bounds.eventgroup(1).lower = [bounds.eventgroup(1).lower,0];
            bounds.eventgroup(1).upper = [bounds.eventgroup(1).upper,0];
        end
        
        if auxdata.settings.default_capital_tax
            bounds.eventgroup(1).lower = [bounds.eventgroup(1).lower, 0,0];
            bounds.eventgroup(1).upper = [bounds.eventgroup(1).upper, 0,0];
        end
        
        if auxdata.settings.equipment_equal_robot_tax
            bounds.eventgroup(1).lower = [bounds.eventgroup(1).lower, 0];
            bounds.eventgroup(1).upper = [bounds.eventgroup(1).upper, 0];
        end
        
        
        
        %-------------------------------------------------------------------------%
        %---------------------- Impose Guesses for Solution ----------------------%
        %-------------------------------------------------------------------------%
        
        % intGuess                            =      -guess.fval;
        intGuess                            =      1./auxdata.settings.inpscale.welfare;
        
        if auxdata.settings.consistency_conditions
            
            guess.phase.integral = [...
                intGuess,...
                intGuess,...
                intGuess,...
                intGuess,...
                0,...
                par_dist.Y_M_0*par0.L_M,...
                par_dist.Y_R_0*par0.L_R,...
                par_dist.Y_C_0*par0.L_C,...
                intGuess,...
                intGuess,...
                intGuess];
        else
            
            guess.phase.integral = [...
                intGuess,...
                intGuess,...
                intGuess,...
                intGuess,...
                0];
            
        end
        
        uGuess                              =      [u0; uf];
        
        lGuess                              =      [l0; lf];
        
        guess.phase.time                    =      [w0; wf];
        guess.phase.state                   =      uGuess;
        guess.phase.control                 =      lGuess;
        
        
        %-------------------------------------------------------------------------%
        %----------Provide Mesh Refinement Method and Initial Mesh ---------------%
        %-------------------------------------------------------------------------%
        mesh.method                         =      'hp-LiuRao-Legendre';
        mesh.tolerance                      =      1e-6;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        mesh.maxiterations                   =     1;
        mesh.colpointsmin                   =      4;
        mesh.colpointsmax                   =      10;
        mesh.phase(1).colpoints = [4, 4, 5, 4, 5, 4, 4, 5, 7, 10];
        
        
        %-------------------------------------------------------------------------%
        %------------- Assemble Information into Problem Structure ---------------%
        %-------------------------------------------------------------------------%
        setup.name                          =      'GPOPS-Problem';
        setup.functions.continuous          =      @GPOPSContinuous;
        setup.functions.endpoint            =      @GPOPSEndpoint;
        setup.bounds                        =      bounds;
        setup.guess                         =      guess;
        setup.mesh                          =      mesh;
        setup.nlp.solver                    =      'ipopt';
        setup.nlp.ipoptoptions.tolerance    =      1e-7;
        setup.nlp.ipoptoptions.linear_solver =     'ma57';
        setup.derivatives.supplier          =      'sparseCD';
        setup.derivatives.derivativelevel   =      'second';
        setup.derivatives.dependencies      =      'full';
        setup.method                        =      'RPM-Differentiation';
        setup.displaylevel                  =      2;
        setup.scales.method                 =      'none';
        % setup.scales.method                 =      'automatic-bounds';
        setup.nlp.ipoptoptions.maxiterations = 200;
        
        auxdata.mesh = setup.mesh;
        
        setup.auxdata                       =      auxdata;
        
        
        % clear results file
        vars_file = [dir_out,auxdata.settings.version,'_vars_',suffix,'.csv'];
        vectors_file = [dir_out,auxdata.settings.version,'_vectors_',suffix,'.csv'];
        if exist(vars_file, 'file')==2
            delete(vars_file);
        end
        
        if exist(vectors_file, 'file')==2
            delete(vectors_file);
        end
        
        %% Handling of parameter bounds
        
        % tolerance for static parameter bounds in terms of guess
        tolParPctLower = 0.5;
        tolParPctUpper = 2;
        
        index_parameter = 1;
           
        eta_R = setup.auxdata.gpops_fun.prod_fun.eta_R;
        kappa_R = setup.auxdata.gpops_fun.prod_fun.kappa_R;
        eta_R_range = linspace(eta_R, 0.75 * eta_R, num_points);
        
        if kappa_R_increase
            kappa_R_range = linspace(kappa_R, 1.25 * kappa_R, num_points);
        else
            kappa_R_range = linspace(kappa_R, 0.75 * kappa_R, num_points);            
        end
       
        
        [ETA_R, KAPPA_R] = meshgrid(eta_R_range, kappa_R_range);
        
        % change the order of every second row of KAPPA_R, such
        % that there are no larger jumps in parameters when
        % linearly walking through them
        
        KAPPA_R(:,[2:2:end]) = flipud(KAPPA_R(:,[2:2:end]));
        
        
        while index_parameter <= numel(KAPPA_R)
            
            %% Load setup after first iteration
            if index_parameter > 1
                setup = load_setup_fun(['output_', suffix, '.mat'])
                eta_R_prev = setup.auxdata.gpops_fun.prod_fun.eta_R;
                kappa_R_prev = setup.auxdata.gpops_fun.prod_fun.kappa_R;
            end
            
            eta_R = ETA_R(index_parameter);
            kappa_R = KAPPA_R(index_parameter);

            % set bounds
            setup = set_bounds(tolParPctLower, tolParPctUpper, setup)
            bounds_violated = 1;
            index_parameter_tmp = index_parameter;
            while bounds_violated
                
                %% Call GPOPS
                nlpinfo = -1;
                while nlpinfo ~= 0
                        setup.auxdata.gpops_fun.prod_fun.eta_R = eta_R;
                        setup.auxdata.gpops_fun.prod_fun.kappa_R = kappa_R;

                        fprintf('\n=============================================================\n')
                        fprintf('running with eta_R = %f, kappa_R = %f\n', eta_R, kappa_R)
                        fprintf('at %s\n', datestr(now))
                        fprintf('=============================================================\n')
                        
                    tic
                    output = gpops2(setup);
                    toc
                    
                    nlpinfo = output.result.nlpinfo;
                    
                    if nlpinfo ~= 0

                        fprintf('\n=============================================================\n')
                        fprintf('failed to solve, reducing step size\n')
                        fprintf('=============================================================\n')
                        
                        eta_R_step = eta_R - eta_R_prev;
                        kappa_R_step = kappa_R - kappa_R_prev;
                        
                        % cut step wrt. to previous value in half
                        eta_R = eta_R_prev + eta_R_step/2;
                        kappa_R = kappa_R_prev + kappa_R_step/2;
                        
                        % reduce index by one, such that once solved, we
                        % continue with the index that caused trouble
                        index_parameter = index_parameter_tmp - 1;
                    end
                
                end
                %% Post-processing
                [stats, vars, vectors] = GPOPSRobustnessComputeVars(output, setup);
                stats
                
                bounds_violated = check_bounds(stats)
                
                if bounds_violated
                    fprintf('before updating bounds\n')
                    setup.bounds.parameter.lower
                    setup.bounds.parameter.upper
                    
                    setup = output.result.setup;
                    
                    setup.bounds.parameter.lower = max(0, 0.5.*setup.bounds.parameter.lower);
                    setup.bounds.parameter.upper = 2.* ...
                        setup.bounds.parameter.upper;
                    
                    fprintf('after updating bounds\n')
                    setup.bounds.parameter.lower
                    setup.bounds.parameter.upper
                end
            end
            
            save(['output_', suffix, '.mat'],'output')
            save_results(vars, vectors);
            index_parameter = index_parameter + 1;
        end
        
        
    end
    diary off

    tools.struct2csv(vars,'solution/vars.csv')
    tools.struct2csv(vectors,'solution/vectors.csv')


    function bounds_violated = check_bounds(guess)
        bounds_violated = false;
        names = fieldnames(guess);
        checks = strncmpi(names,'check',5);
        check_names = names(checks);
        
        n = length(check_names);
        index_bound_check = 1;
        while ~bounds_violated && index_bound_check <= n
            if guess.(cell2mat(check_names(index_bound_check))) < 1e-2
                fprintf('bound for %s violated', ...
                        cell2mat(check_names(index_bound_check)))
                bounds_violated = true;
            end
            index_bound_check = index_bound_check + 1;
        end
    end

    function setup = load_setup_fun(file)
        fprintf('loading setup\n')
        clear('setup')
        load(file)
        setup = output.result.setup;
        auxdata = setup.auxdata;
        
        par = output.result.solution.parameter;
        
        % adjust inpscale
        inpscale = setup.auxdata.settings.inpscale;
        
        inpscale.L_M = par(1) * inpscale.L_M;
        inpscale.L_R = par(2) * inpscale.L_R;
        inpscale.L_C = par(3) * inpscale.L_C;
        inpscale.K_B_M = par(4) * inpscale.K_B_M;
        inpscale.K_B_R = par(5) * inpscale.K_B_R;
        inpscale.K_B_C = par(6) * inpscale.K_B_C;
        inpscale.K_E = par(7) * inpscale.K_E;
        inpscale.K_S = par(8) * inpscale.K_S;
        inpscale.transfer = par(9) * inpscale.transfer;
        
        setup.auxdata.settings.inpscale = inpscale;
        
        % adjust guess and bounds
        setup.guess.parameter = ones(1,9);
        
        setup.bounds.parameter.lower(:) = 0.5;
        setup.bounds.parameter.upper(:) = 2;
        
        setup.bounds.parameter.upper(4) = 100;
        setup.bounds.parameter.upper(5) = 100;
        setup.bounds.parameter.upper(6) = 100;
    end

    function setup = set_bounds(tolParPctLower, tolParPctUpper, setup)
        
        guess.parameter = ones(1,9);
        
        bounds.parameter.lower = max(guess.parameter - tolParPctLower.*guess.parameter,0);
        bounds.parameter.upper = guess.parameter + tolParPctUpper.*guess.parameter;
        
        bounds.parameter.upper(5) = 100;
        bounds.parameter.lower(9) = 1e-3;
        bounds.parameter.upper(9) = 100;
        
        setup.bounds.parameter = bounds.parameter;
        setup.guess.parameter = guess.parameter;
        
    end

    function save_results(vars, vectors)
        try
            % load results file, provided it exists
            results_vars = struct2table(tools.csv2struct(file_vars));
            results_vectors = struct2table(tools.csv2struct(file_vectors));
        end
        
        results_vars = [results_vars; struct2table(vars)];
        results_vectors = [results_vectors; struct2table(vectors)]; ...
            
        results_vars_struct = table2struct(results_vars);
        results_vectors_struct = table2struct(results_vectors);
        
        tools.struct2csv(results_vars_struct, file_vars)
        tools.struct2csv(results_vectors_struct, file_vectors)
    end

end % sub-functions defined inside main function to be able to use
    % arguments that are not explicitly passed